Input and Output Paths

if(exists("snakemake")) {
  input_list <- snakemake@input
  output_list <- snakemake@output
} else {
  # inputs:
  input_list <- list(
    oregon_rivers = "inputs/Rivers_OR/Rivers_OR.shp"
  )
  # outputs:
  output_list <- list(
    texmap = "tex/images/map-crop.pdf"
  )
}

# add the other entries based on the final texmap
output_list$map = "results/map_of_samples/map.pdf"
output_list$cropmap = "results/map_of_samples/map-crop.pdf"


# we create the necessary output directories like this:
dump <- lapply(output_list, function(x)
  dir.create(dirname(x), recursive = TRUE, showWarnings = FALSE)
)
library(terra)
library(sf)
library(tidyverse)
library(ggspatial)
#library(ggsn)
library(cowplot)
library(plotly)
library(maps) # needed for map_data()
library(mapproj)  # making it explicit for renv

Download spatial data

Raster

if(!file.exists("geo-spatial/HYP_HR_SR_OB_DR/HYP_HR_SR_OB_DR.tif")) {
  dir.create("geo-spatial", showWarnings = FALSE)
  download.file(
    url = "https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/raster/HYP_HR_SR_OB_DR.zip",
    destfile = "geo-spatial/HYP_HR_SR_OB_DR.zip"
  )
  
  unzip(
    zipfile = "geo-spatial/HYP_HR_SR_OB_DR.zip", 
    exdir = "geo-spatial/HYP_HR_SR_OB_DR"
  )
  
  file.remove("geo-spatial/HYP_HR_SR_OB_DR.zip")
}

State lines

if(!file.exists("geo-spatial/ne_10m_admin_1_states_provinces_lines/ne_10m_admin_1_states_provinces_lines.shp")) {
  dir.create("geo-spatial", showWarnings = FALSE)
  download.file(
    url = "https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/cultural/ne_10m_admin_1_states_provinces_lines.zip",
    destfile = "geo-spatial/ne_10m_admin_1_states_provinces_lines.zip"
  )
  
  unzip(
    zipfile = "geo-spatial/ne_10m_admin_1_states_provinces_lines.zip", 
    exdir = "geo-spatial/ne_10m_admin_1_states_provinces_lines"
  )
  
  file.remove("geo-spatial/ne_10m_admin_1_states_provinces_lines.zip")
}

Coastline

if(!file.exists("geo-spatial/ne_10m_coastline/ne_10m_coastline.shp")) {
  dir.create("geo-spatial", showWarnings = FALSE)
  download.file(
    url = "https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/physical/ne_10m_coastline.zip",
    destfile = "geo-spatial/ne_10m_coastline.zip"
  )
  
  unzip(
    zipfile = "geo-spatial/ne_10m_coastline.zip", 
    exdir = "geo-spatial/ne_10m_coastline"
  )
  
  file.remove("geo-spatial/ne_10m_coastline.zip")
}

California Major Rivers and Creeks

This is from California Nat Res Agency. Find it at: https://data.cnra.ca.gov/dataset/national-hydrography-dataset-nhd

# note.  Once when I did this, download.file did not see capable of downloading this file on my Mac
# (the download is super sloooooow), and that time I 
# ended up downloading it with Chrome.  But it seems those connection problems
# had been fixed when I did this again.
if(!file.exists("geo-spatial/NHD_Major_Rivers_and_Creeks/Major_Rivers_and_Creeks.shp")) {
  dir.create("geo-spatial", showWarnings = FALSE)
  download.file(
    url = "https://data.cnra.ca.gov/dataset/511528b2-f7d3-4d86-8902-cc9befeeeed5/resource/7d1e7e44-81b1-43fe-95f6-1862eea6ac24/download/nhd_major_rivers_and_creeks.zip",
    destfile = "geo-spatial/nhd_major_rivers_and_creeks.zip"
  )
  
  unzip(
    zipfile = "geo-spatial/nhd_major_rivers_and_creeks.zip", 
    exdir = "geo-spatial"
  )
  
  file.remove("geo-spatial/nhd_major_rivers_and_creeks.zip")
}

Read in the spatial data

nat.earth <- rast("geo-spatial/HYP_HR_SR_OB_DR/HYP_HR_SR_OB_DR.tif") #
state_prov <- st_read("geo-spatial/ne_10m_admin_1_states_provinces_lines/ne_10m_admin_1_states_provinces_lines.shp")
## Reading layer `ne_10m_admin_1_states_provinces_lines' from data source 
##   `/Users/eriq/Documents/git-repos/california-chinook-microhaps/geo-spatial/ne_10m_admin_1_states_provinces_lines/ne_10m_admin_1_states_provinces_lines.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 10179 features and 57 fields (with 1 geometry empty)
## Geometry type: MULTILINESTRING
## Dimension:     XY
## Bounding box:  xmin: -178.1371 ymin: -49.25087 xmax: 178.4486 ymax: 81.12853
## Geodetic CRS:  WGS 84
coastline <- st_read("geo-spatial/ne_10m_coastline/ne_10m_coastline.shp")
## Reading layer `ne_10m_coastline' from data source 
##   `/Users/eriq/Documents/git-repos/california-chinook-microhaps/geo-spatial/ne_10m_coastline/ne_10m_coastline.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 4133 features and 3 fields
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: -180 ymin: -85.22194 xmax: 180 ymax: 83.6341
## Geodetic CRS:  WGS 84
all_rivers <- st_read("geo-spatial/NHD_Major_Rivers_and_Creeks/Major_Rivers_and_Creeks.shp") %>%
  st_zm() %>%
  st_transform(., st_crs(state_prov))
## Reading layer `Major_Rivers_and_Creeks' from data source 
##   `/Users/eriq/Documents/git-repos/california-chinook-microhaps/geo-spatial/NHD_Major_Rivers_and_Creeks/Major_Rivers_and_Creeks.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 220702 features and 17 fields
## Geometry type: LINESTRING
## Dimension:     XYZM
## Bounding box:  xmin: -124.4252 ymin: 32.50998 xmax: -113.9205 ymax: 42.6784
## z_range:       zmin: -2e-04 zmax: 1476.791
## m_range:       mmin: -1.797693e+308 mmax: 100
## Geodetic CRS:  NAD83 + NAVD88 height
or_rivers <- st_read(input_list$oregon_rivers)
## Reading layer `Rivers_OR' from data source 
##   `/Users/eriq/Documents/git-repos/california-chinook-microhaps/inputs/Rivers_OR/Rivers_OR.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 83 features and 8 fields
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: -124.4226 ymin: 38.65888 xmax: -108.1906 ymax: 49.27932
## Geodetic CRS:  WGS 84

Crop the basemap, state lines coast lines

# important to put them in this order and named like this
domain <- c(
  xmin = -127.9,
  xmax = -117.4,
  ymin = 36,
  ymax = 43.25
)
sf_use_s2(FALSE)
nat_crop <- crop(nat.earth, y = ext(domain))
state_subset <- st_crop(state_prov, domain)
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
coastline_cropped <- st_crop(coastline, domain)
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries

Subset the rivers

We plot this with plotly so we can see the names of the rivers easily, and whittle it down to the ones that we need.

# I don't know the gnis_id for mill creek, and there are a lot of Mill Creeks,
# so lets get it:
mill_area <- c(
  xmin = -121.4,
  xmax = -122.3,
  ymin = 40,
  ymax = 40.3
)
mill_candi <- all_rivers %>%
  filter(str_detect(gnis_name, "Mill Creek")) %>%
  st_crop(mill_area)
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
# same with Blue Creek
blue_area <- c(
  xmin = -124,
  xmax = -123.5,
  ymin = 41.2,
  ymax = 41.8
)

blue_candi <- all_rivers %>%
  filter(str_detect(gnis_name, "Blue Creek")) %>%
  st_crop(blue_area)
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
# rogue river:
rogue <- or_rivers %>%
  filter(name == "Rogue") %>%
  mutate(gnis_name = "Rogue River")


our_rivers <- all_rivers %>%
  filter(gnis_name %in% c(
    "Sacramento River",
    "San Joaquin River",
    "Feather River",
    "Russian River",
    "Trinity River",
    "Klamath River",
    "Eel River",
    "Smith River"
    #"Deer Creek",
    #"Battle Creek"
  ) | 
    gnis_id %in% c(
      "00220293", "00237425", "00263498", "00266522",  # these are butte creek
      "00233775", "01655075"    # these are deer
     # "00218740", "00229640", "00234966"   # these are battle creek
    ) 
  ) %>%
  bind_rows(mill_candi) %>%
  bind_rows(blue_candi) %>%
  bind_rows(rogue)


# get candidate positions for notations
labels <- read_tsv("inputs/map-notations.tsv", comment = "#")


g <- ggplot() + 
  geom_sf(data = our_rivers, aes(colour = gnis_name)) +
  geom_sf(data = coastline_cropped) +
  geom_point(data = labels, aes(x = label_long, y = label_lat), colour = "red") +
  geom_point(data = labels, aes(x = tip_long, y = tip_lat), colour = "blue", size = 0.2)

ggplotly(g)

That looks good. So let us proceed.

Make the base map

base_map <- ggplot() +
  ggspatial::layer_spatial(nat_crop) +
  geom_sf(data = state_subset, color = "gray30", fill = NA) +
  geom_sf(data = coastline_cropped, color = "gray30", fill = NA, linewidth = 0.15) +
  geom_sf(data = our_rivers, colour = "blue", linejoin = "round", lineend = "round")

Now add stuff to that:

bit <- 0.0
kick <- 0.4
source("R/colors.R")

labels_1line <- labels %>%
  group_by(map_text) %>%
  filter(rank == max(rank)) %>%
  ungroup()

# define the shapes
point_shapes <- c(
  `Fall run` = 21,
  `Spring run` = 22,
  `Winter run` = 23,
  `Late-fall run` = 24
)

mm <- base_map + 
  geom_segment(data = labels, aes(x = tip_long, y = tip_lat, xend = label_long, yend = label_lat), colour = "black", linewidth = 0.4) +
  geom_point(data = labels, aes(x = label_long + bit + (1 - hjust * 2) * kick * (rank - 1), y = label_lat, fill = run_timing, shape = run_timing), size = 3.5) +
  geom_label(
    data = labels_1line,
    aes(
      x = label_long + bit + (1 - hjust * 2) * kick * ( 0.65 + rank - 1),
      y = label_lat, 
      label = map_text,
      hjust = hjust
    ),
    size = 2.2,
    label.padding = unit(0.09, "lines"),
    lineheight = 0.85
  ) +
  scale_fill_manual(values = run_time_colors) +
  scale_shape_manual(values = point_shapes) +
  theme_bw() +
  theme(
    panel.border = element_rect(colour = "black", linewidth = 1),
    axis.text.x = element_text(size = 8, family = "serif", angle = 35, hjust = 1),
    axis.text.y = element_text(size = 8, family = "serif"),
    axis.title.y = element_text(family = "serif", size = 10),
    axis.title.x = element_text(family = "serif", vjust = 2, size = 10),
    plot.margin = margin(0, 0.1, 0, 0.15, "cm"),
    legend.position = "none"
  ) +
  xlab("Longitude") +
  ylab("Latitude") +
  coord_sf(
    expand = FALSE,
  ) +
  guides(fill = guide_legend(title = "Run Timing:", nrow = 2)) +
  ggspatial::annotation_north_arrow(location = "br", height = unit(10, "mm"), width = unit(5, "mm"), style = north_arrow_fancy_orienteering()) +
  ggspatial::annotation_scale(location = "br", height = unit(0.1, "cm"))

Now, work on the world-scale map with the inset:

wrld <- map_data("world")
states <- map_data("state")
domain_df <- tibble(point = 1:length(domain), long = rep(domain[1:2], each = 2), lat = c(domain[3:4], rev(domain[3:4])))

inset_world <- ggplot() +
  geom_polygon(data = wrld, aes(x = long, y = lat, group = group), colour = "black", fill = "gray90", linewidth = 0.1) +
  geom_path(data = states, aes(x = long, y = lat, group = group), colour = "black", linewidth = 0.1) +
  geom_polygon(data = domain_df, mapping = aes(x = long, y = lat), colour = "red", fill = "red", alpha = 0.3, linewidth = 0.1) +
  coord_map("ortho", orientation = c(30, -114, 0), xlim = c(-150, -75), ylim = c(28, 60)) +
  theme_bw() +
  labs(x = NULL, y = NULL) +
  theme(
    axis.title = element_blank(),
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    plot.margin = unit(c(0, 0, -1, -1), "mm")
  )
inset_world

Use cowplot to make the inset appear on the bigger map

final_map <- ggdraw() +
  draw_plot(mm) +
  draw_plot(inset_world, x = 0.215, y = 0.139, width = 0.20, height = 0.17)

ggsave(final_map, filename = output_list$map, width = 5, height = 3.5)

# let's crop that down for use, too
CALL <- paste("pdfcrop ", output_list$map, collapse = " ")
system(CALL)

file.copy(from = output_list$cropmap, to = output_list$texmap, overwrite = TRUE)
## [1] TRUE

Make colored balls for inclusion into the caption

We can do this and just kick them directly into the tex/images folder.

rc_col_tib <- enframe(run_time_colors) %>%
  mutate(
    nospace = str_replace_all(name, " ", "-"),
    file = str_c("tex/images/", nospace, "-ball.pdf")
  )

for(r in 1:nrow(rc_col_tib)) {
  tmp <- rc_col_tib[r,]
  g <- ggplot(tmp, aes(x = 1, y = 1)) +
    geom_point(shape = point_shapes[tmp$name], fill = tmp$value, size = 3.5) +
    theme_void()
  ggsave(g, filename = tmp$file)
  
  CALL <- str_c("pdfcrop ", tmp$file)
  system(CALL)
  file.remove(tmp$file)
}

Session Info

sessioninfo::session_info()
## ─ Session info ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
##  setting  value
##  version  R version 4.4.3 (2025-02-28)
##  os       macOS Monterey 12.7.5
##  system   aarch64, darwin20
##  ui       X11
##  language (EN)
##  collate  en_US.UTF-8
##  ctype    en_US.UTF-8
##  tz       America/Denver
##  date     2025-03-27
##  pandoc   3.6.4 @ /Users/eriq/Documents/git-repos/california-chinook-microhaps/.snakemake/conda/def120e971fb2c87a8c042b4c2c09bdb_/bin/ (via rmarkdown)
##  quarto   1.4.550 @ /usr/local/bin/quarto
## 
## ─ Packages ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
##  ! package     * version date (UTC) lib source
##  P bit           4.6.0   2025-03-06 [?] CRAN (R 4.4.1)
##  P bit64         4.6.0-1 2025-01-16 [?] CRAN (R 4.4.1)
##  P bslib         0.9.0   2025-01-30 [?] CRAN (R 4.4.1)
##  P cachem        1.1.0   2024-05-16 [?] CRAN (R 4.4.1)
##  P class         7.3-23  2025-01-01 [?] CRAN (R 4.4.3)
##  P classInt      0.4-11  2025-01-08 [?] CRAN (R 4.4.1)
##  P cli           3.6.4   2025-02-13 [?] CRAN (R 4.4.1)
##  P codetools     0.2-20  2024-03-31 [?] CRAN (R 4.4.3)
##  P colorspace    2.1-1   2024-07-26 [?] CRAN (R 4.4.1)
##  P cowplot     * 1.1.3   2024-01-22 [?] CRAN (R 4.4.0)
##  P crayon        1.5.3   2024-06-20 [?] CRAN (R 4.4.1)
##  P crosstalk     1.2.1   2023-11-23 [?] CRAN (R 4.4.0)
##  P data.table    1.17.0  2025-02-22 [?] CRAN (R 4.4.1)
##  P DBI           1.2.3   2024-06-02 [?] CRAN (R 4.4.1)
##  P digest        0.6.37  2024-08-19 [?] CRAN (R 4.4.1)
##  P dplyr       * 1.1.4   2023-11-17 [?] CRAN (R 4.4.0)
##  P e1071         1.7-16  2024-09-16 [?] CRAN (R 4.4.1)
##  P evaluate      1.0.3   2025-01-10 [?] CRAN (R 4.4.1)
##  P farver        2.1.2   2024-05-13 [?] CRAN (R 4.4.1)
##  P fastmap       1.2.0   2024-05-15 [?] CRAN (R 4.4.1)
##  P forcats     * 1.0.0   2023-01-29 [?] CRAN (R 4.4.0)
##  P generics      0.1.3   2022-07-05 [?] CRAN (R 4.4.1)
##  P ggplot2     * 3.5.1   2024-04-23 [?] CRAN (R 4.4.0)
##  P ggspatial   * 1.1.9   2023-08-17 [?] CRAN (R 4.4.0)
##  P glue          1.8.0   2024-09-30 [?] CRAN (R 4.4.1)
##  P gtable        0.3.6   2024-10-25 [?] CRAN (R 4.4.1)
##  P hms           1.1.3   2023-03-21 [?] CRAN (R 4.4.0)
##  P htmltools     0.5.8.1 2024-04-04 [?] CRAN (R 4.4.1)
##  P htmlwidgets   1.6.4   2023-12-06 [?] CRAN (R 4.4.0)
##  P httr          1.4.7   2023-08-15 [?] CRAN (R 4.4.0)
##  P jquerylib     0.1.4   2021-04-26 [?] CRAN (R 4.4.0)
##  P jsonlite      1.9.1   2025-03-03 [?] CRAN (R 4.4.1)
##  P KernSmooth    2.23-26 2025-01-01 [?] CRAN (R 4.4.3)
##  P knitr         1.49    2024-11-08 [?] CRAN (R 4.4.1)
##  P labeling      0.4.3   2023-08-29 [?] CRAN (R 4.4.1)
##  P lazyeval      0.2.2   2019-03-15 [?] CRAN (R 4.4.1)
##  P lifecycle     1.0.4   2023-11-07 [?] CRAN (R 4.4.1)
##  P lubridate   * 1.9.4   2024-12-08 [?] CRAN (R 4.4.1)
##  P magrittr      2.0.3   2022-03-30 [?] CRAN (R 4.4.1)
##  P mapproj     * 1.2.11  2023-01-12 [?] RSPM
##  P maps        * 3.4.2.1 2024-11-10 [?] RSPM
##  P munsell       0.5.1   2024-04-01 [?] CRAN (R 4.4.1)
##  P pillar        1.10.1  2025-01-07 [?] CRAN (R 4.4.1)
##  P pkgconfig     2.0.3   2019-09-22 [?] CRAN (R 4.4.1)
##  P plotly      * 4.10.4  2024-01-13 [?] CRAN (R 4.4.0)
##  P proxy         0.4-27  2022-06-09 [?] CRAN (R 4.4.1)
##  P purrr       * 1.0.4   2025-02-05 [?] CRAN (R 4.4.1)
##  P R6            2.6.1   2025-02-15 [?] CRAN (R 4.4.1)
##  P ragg          1.3.3   2024-09-11 [?] CRAN (R 4.4.1)
##  P Rcpp          1.0.14  2025-01-12 [?] CRAN (R 4.4.1)
##  P readr       * 2.1.5   2024-01-10 [?] CRAN (R 4.4.0)
##    renv          1.1.4   2025-03-20 [1] CRAN (R 4.4.1)
##  P rlang         1.1.5   2025-01-17 [?] CRAN (R 4.4.1)
##  P rmarkdown     2.29    2024-11-04 [?] CRAN (R 4.4.1)
##  P sass          0.4.9   2024-03-15 [?] CRAN (R 4.4.0)
##  P scales        1.3.0   2023-11-28 [?] CRAN (R 4.4.0)
##  P sessioninfo   1.2.3   2025-02-05 [?] CRAN (R 4.4.1)
##  P sf          * 1.0-20  2025-03-24 [?] RSPM (R 4.4.0)
##  P stringi       1.8.4   2024-05-06 [?] CRAN (R 4.4.1)
##  P stringr     * 1.5.1   2023-11-14 [?] CRAN (R 4.4.0)
##  P systemfonts   1.2.1   2025-01-20 [?] CRAN (R 4.4.1)
##  P terra       * 1.8-29  2025-02-26 [?] CRAN (R 4.4.1)
##  P textshaping   1.0.0   2025-01-20 [?] CRAN (R 4.4.1)
##  P tibble      * 3.2.1   2023-03-20 [?] CRAN (R 4.4.0)
##  P tidyr       * 1.3.1   2024-01-24 [?] CRAN (R 4.4.1)
##  P tidyselect    1.2.1   2024-03-11 [?] CRAN (R 4.4.0)
##  P tidyverse   * 2.0.0   2023-02-22 [?] CRAN (R 4.4.0)
##  P timechange    0.3.0   2024-01-18 [?] CRAN (R 4.4.1)
##  P tzdb          0.4.0   2023-05-12 [?] CRAN (R 4.4.0)
##  P units         0.8-7   2025-03-11 [?] CRAN (R 4.4.1)
##  P vctrs         0.6.5   2023-12-01 [?] CRAN (R 4.4.0)
##  P viridisLite   0.4.2   2023-05-02 [?] CRAN (R 4.4.1)
##  P vroom         1.6.5   2023-12-05 [?] CRAN (R 4.4.0)
##  P withr         3.0.2   2024-10-28 [?] CRAN (R 4.4.1)
##  P xfun          0.51    2025-02-19 [?] CRAN (R 4.4.1)
##  P yaml          2.3.10  2024-07-26 [?] CRAN (R 4.4.1)
## 
##  [1] /Users/eriq/Documents/git-repos/california-chinook-microhaps/renv/library/macos/R-4.4/aarch64-apple-darwin20
##  [2] /Users/eriq/Library/Caches/org.R-project.R/R/renv/sandbox/macos/R-4.4/aarch64-apple-darwin20/f7156815
## 
##  * ── Packages attached to the search path.
##  P ── Loaded and on-disk path mismatch.
## 
## ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────

Running Time

Running the code and rendering this notebook required approximately this much time:

Sys.time() - start_time
## Time difference of 16.5015 secs